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Abstract 

We introduce fractional Brownian motion processes (fBm) as an alternative model 
for the turbulent index of refraction. These processes allow to reconstruct most 
of the refractive index properties, but they are not differentiable. We overcome 
the apparent impossibility of their use within the Ray Optics approximation in- 
troducing a Stochastic Calculus. Afterwards, we successfully provide a solution for 
the stochastic ray-equation; moreover, its implications in the statistical analysis 
of experimental data is discussed. In particular, we analyze the dependence of 
the averaged solution against the characteristic variables of a simple propagation 
problem. 
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I. INTRODUCTION 



Diverse experimental techniques have been devoted to the study of the optical properties 
of the turbulent atmosphere. Most of these techniques are based on the analysis of the 
output of laser beams making their way through it. But also, controlled experiences had 
been developed for the laboratory, such as the experiments performed by Consortini et al^2£. 
These experiences apply Geometric Optics to interpret the acquired data. This analysis has 
its theoretical grounds on the precursor paper by P. Beckmar4, who was able to find a simple 
relationship between the variance of the turbulent refractive index /x(r) — being homogeneous 
and isotropic — and the variance of the laser beam wandering over a screen. This derivation, 
however, is based upon the strong assumption that the stochastic process is smooth enough 
to allow continuity in its derivatives. 

This assumption has serious problems when checked against Tatarski's foundational work 
about lightwave propagation through turbulent atmosphere^. A short revision of the covari- 
ance of the turbulent index of refraction found by him shows that the stochastic process 
associated is nowhere differentiable 6 . 

However, the model introduced by Beckmann is just one of the many used during the 
last 40 years to solve imaging problems through the atmosphere, with relative success. 
In fact, it is clear that the stochastic properties of the refractive index in Atmospheric 
Optics have never been fully understood nor explained. For instance, most works treated it 
like a Gaussian process^, some others suggested stationary increments, while other works 
have proposed the use of non-Gaussian statistics. Ishimaru's book^ presents an extensive 
description of these works. 

Simultaneously, over the last decade an intense debate in Fluid Dynamics has been carried 
out about whether or not passive scalar fields, among which is the turbulent refractive 
index, behave like the velocity field. That is, under which circumstances they inherit the 
stochastic properties of the turbulent velocity u(r), which is a Gaussian process that follows 
the Kolmogorov refined similarity hypotheses inside the inertial range, Iq -C ||r|| <C Lq, 



where e r is the dissipated energy, A n is a constant depending on n; also, lo (inner length) and 
L (outer length) are constants that determine the inertial range. They can be estimated the- 
oretically. When intermittence effects are noticeable the dissipated energy modifies slightly 




(1) 
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the right-hand side of this equation, changing the power over the correlation distance 

n n 

3 - 3 + Cn 

here ( n is called multi-fractal exponent. 

It was show n 10 ! 11 that passive scalars fields are nearly gaussian as far as Iq <C L q , and are 
unsatisfactory modeled either by log-normal distributions, or the so called Frisch's /3-model. 
Effectively, in the inertial range, these fields obey a law resembling the Kolmogorov's law 
for the velocity fields. Moreover, they do not present the multi-fractal property due to 
intermittence whenever an isotropic velocity field is present^, i.e. ( n = 0. Hence, the 
behavior of the turbulent refractive index predicted by Tatarski has been confirmed and 
extended. 

In the meantime, Stolovitzky and Sreenivasan^ successfully obtained Eq. (£[} modeling 
the turbulent velocity field as a fractional Brownian motion (fBm). But this model failed 
to replicate the intermittent property. It must be stressed that Kolmogorov refined simi- 
larity hypotheses also implies that the velocity field is independent of the dissipated energy 
probability distribution. 

According to what we have pointed out here, the fBm processes seem to be a good al- 
ternative model for the turbulent refractive index. First, they are gaussian, and second, 
they let us test the Structure Function's power factor. Also, they are continuos but nowhere 
differentiable, as it results from the application of the Kolmogorov hypotheses to the refrac- 
tive index. In few words, the fractional Brownian motion model for the turbulent index of 
refraction describes closer the turbulent properties of passive scalar fields than Beckmann's 
proposed model. However, a mathematical difficulty is introduced since we have lost dif- 
ferentiability in the usual sense, and as a consequence we cannot appeal to the variational 
methods used in Geometric Optics. Our task here will be to show how using Stochastic 
Calculus techniques (for a first introduction we refer the reader to 0ksenda l 14il fr ) this sit- 
uation can be overridden providing an explicit solution to the ray propagation through air 
in turbulent motion. Besides the intrinsic complexity of these tools, our model is meant to 
provide a bridge between the stochastic processes and experimental data analysis. Also, we 
will gain knowledge about problems handling (geodesic) stochastic equations in Optics. 
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II. STOCHASTIC DIFFERENTIAL EQUATIONS IN GEOMETRIC OPTICS 



As it is well know n 16 ' 17 , the Lagragian inherited from the application of the Fermat's 
Extremal Principle in Geometric Optics is singular. In the Appendix [X] we show the ray- 
equations associated to this Lagragian result to be, Eqs. (|A5J) : 

d ( 1 dq \ X^ 2 f^.f_\\ j xi i.-_„ || 1 1 2 \2 2, 



, . "V„n ((?(r)), and the constraint equation \\q\\ = X n (g), 

dr \Xdr J 2 

where q(r) : R — > M 3 is the ray-light trajectory with parameter r, n is the refractive index, 
while the smooth function A : M. 3 — > IR can be freely choosen. The election of this function 
fixes the parametrization. 

Since the equations above are nonlinear, we will begin our study by linearizing them. 
Moreover, at the same time we will set the parameter r. Let n be the refractive index of 
the medium and no its average, we write 

n 2 (q) =n 2 + ae(q), (2) 

e(q) represents a perturbation field with its intensity measured by a. We also assume that 
it contains all the inhomogeneities of the media, so when a = the index is constant. Now 
we express the solution to Eq. ()A5|1 in a power series on a: 

Q{ r ) = <?o + ocqi + a 2 q 2 H . 

Although we can also develop a series for the constraint function A, it is far more con- 
venient to set its value beforehand. Let us rewrite the first equation in Eq. (jA5|) as follows 



aX V g e+ — (V 9 A ■ q) — 



dr 2 2 



(3) 



From all the possible parameterization we choose Eq. ()A6j) . so 

1 1 ^ (-l) n e n (q) 



1 1 



n 2 (q) n 2 ' n 2 D n+2 



1 "n 

n=l u 

in short we will write X 2 n := (-l) n t n (q)/n 2 Q n+2 . 

Inserting the former series for A 2 in ©, after some algebraic manipulation we obtain the 
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following family of differential Eqs. 
d 2 q 



dr 2 

dr 2 

d 2 q n 
dr 2 



0. 



2ng 



9 dr J dr 



m=l 



m , 

/c=i ^ 



and when n > 1 : 
dqm—k 



(4) 



dr 



the constraint condition also gives us a constraint equation for each differential equation in 
Eq. ©; 

/ da,* \ ^ 



dg 
dr 
dgi dg 
dr dr 

Edf/fc dg n _fc 
dr dr 

fc=0 



0, 
0. 



for all n > 1. 



We can readily find the zero-order solution of the first equation in (@J). The result is the 
linear relationship: qo(r) = ar + b. Given that the boundary condition to this problem is 



3(0) = 0, 



(5) 



it implies that b = 0. Now we use the constraint condition to obtain ||a|| 2 = 1, so we are 
free to choose the coordinate frame best suited to our purposes. Let this coordinate frame 
be: 

ze z := g = re z , 

this will be our forward direction of propagation. To solve the remaining equations we need 
a bit more than algebra. 

The turbulent refractive index measures the separation between the index of refraction 
and its average; /i(r) := n(r) — no- It is a small quantity, that is why its increments are 
often replaced in the literature by those of the permitivity. This passive scalar field follows 
the Kolmogorov refined hypotheses in the inertial range, and so its Structure Function is 



D £ (r)=E (n 2 (r + r) -n 2 (r')f 



4C 2 f /3 - 2H \\r\\ 2H , 



(6) 



C 2 is the permitivity structure constant and H is some positive constant less than one. If the 
turbulence is isotropic and homogeneous then the Kolmogorov hypotheses sets H = 1/3; so, 
we have introduced the inner scale to correct the departure from this ideal situation. Then, 



according to what we have just said at the introduction and the definition (J2J), we propose 
the following: 

e(r):=B H (l^\\r\\) (7) 

B H is a fractional Brownian motion, and H e (0, 1) is the Hurst parameter. It is a Gaussian 
process with the following properties^: 



E[B H (s)] = 0, 



1 r 

2 



I 'IH i ill i ( 1 2,U 

s\ + \t\ — \s — t\ 



E[B H (s)B H (t)} 
and the self-similarity property 

B H (as) = a H B H (s), for any a, 



(8) 



this last equation implies that both variables have the same probability distribution. Thus 
from Eqs. © and (J7|) we have, 



D e (r) 



4C*l 2 /3 - 2H \\r\\ 2H 



a 2 E 



(B H (lo U 



H/i-l II vi 



B"(r 



..2i-2H\\ r \\2H 
L II 1 II ) 



whenever llrll <C ||r'||, when, 



Estimates for the Structure Constant and the inner length tell us that a ~ 10~ 6 . Therefore, 
in order to examine the stochastic behavior of a wandering beam will be enough to consider 
the first order solution. 

The first order constraint condition reads then 

dql 



dz 



0. 



and together with boundary condition (jHJ) makes the component along the z-axis null all 
over the ray trajectory. Of course, this condition agrees with the corresponding dynamical 
Eq. (J1J. So we are left with a differential equation for the perpendicular displacements to 
the direction of propagation. Finally, we multiply these displacements by a, and write the 
first-order equation as 

d 2 a 

_Q V„ f ( 2 e l + Q). (9) 



We must provide a context to understand the previous equation. That is, a stochastic 
equation is not only determined by the type of process (the fractional Brownian motion in 
our case) attached to it, but also by the integro- differential theory employed to define its 
derivatives. Moreover, there are distinctive stochastic integration methods whether if > 1/2 
or H < 1/2 19 . Here we are going to make use of the Stochastic Calculus exposed in the 
Appendix El so only the H > 1/2 case will be considered. By doing so, either we are 
considering the inert ial- diffusive range, in the following sense 

or the anisotropic scalar situation ( n —>■ l 20 . This situation could be observed in a laboratory 
if an isotropic velocity field can not achieved by the experimental setup. 

Because the turbulent refractive index oscillates around its mean value, it is expected 
that the light wanders around the z-axis over the screen. So the solution we are looking for 
must have expectation zero. This can easily achieved by the formalism we are employing: 
the stochastic integrals (formally known as fractional ltd integrals) defined by the fractional 
white noise and Wick product on fractional Hida spaces have expectation zero. Henceforth, 
from definition (|ZJ) we can calculate the gradient of the index of refraction. We have, using 
the continuity and differential properties for the fractional Brownian motion — Eqs. (|B10Jl 
and ()B12j) — in Sjj and applying the chain rule, the following 

d ,-„„„_, n dB H x l W H (/o 1 "-"^ 



dx l ds 



s= ,-i|| r H ZolMI / || r || 



where W H is the fractional white noise. Remember, once more, this last identity must be 
understood in terms of the formal definition of white noise inside the fractional Hida spaces, 
it has nothing to do with the usual concept of derivative. 

Next, the procedure to interpret Eq. (JHJ) requires to replace all the ordinary products 
containing stochastic variables by Wick products, and we finally write 

'w H (l l \\ze z + Q\\) 



d 2 a 
dz 2 2 l nl 



\ze z + Q\\ 



oQ. (10) 



The fractional white noise is a functional on the real line and its composition with another 
stochastic process has to be defined. Because any analytic function is expressed by a power 
series, as 0ksendal et al. lA suggest, we follow our substitution rule for products and replace 



7 



the powers in the series by Wick's powers — just as we did with the Wick exponential in 
Eq. (jB14j) from the Appendix [B] — whenever a stochastic process is an argument for the 
given function. The representation for the noise in is a series with analytic functions as 
components, Eq. (jBll|) ; thus, we change these components 

(f)(s, z) i k (s) ds -> / (jf{s,Z)l k {s)ds, 
Jr 

Z is some continuous stochastic process with E[Z] := zq ^ 0, and •) is the Wick 
representation of 0(s, •) = H(2H — 1) \s — -\ 2H 2 . We can write \\ze z + Q\\ — z + Q 2 /2z 
because Q ~ 0(a), and then we can evaluate the fractional white noise at z + a 2 Z{uj): 

<P°(s, z + a 2 Z) = H(2H -l)\z + a 2 Z - s | o(2H - 2) = H ( 2H - 1) [{z - s) + a 2 zY (2H ~ 2) , 

we have just took the positive part of the absolute value; it is enough for us examine this 
situation. If E[a 2 Z] = a 2 z then 



[(z -s) + a 2 Z] 



2 --10(2/7-2) 



u 2 l2g -2 , ^ g 2 "(2H-2)---(2H-3-n) 

[( ^ S) + ^° ] + g nl[(z-s) + a ^r 2 - 2 - (Z - Z0) ' 



and all the terms in the series are of order higher or equal to 2 in a. We just need to compare 
the first term against the deterministic coefficient in the white noise series, 

(p(s, t + a 2 z ) - <p(s, t) ~ (t- s) {2H ~ 3) a 2 

and because this happens coordinate to coordinate in the fractional white noise decomposi- 
tion we find 

W H (l^z) W H (l^\\ze z + Q\\) 2 ^ 

- U{a 



z \\ze z + Q\\ 

The first-order Eq. (fTU|) is unaffected by this replacement since they differ in a 2 . We have 
arrived at the linear equation, 



d 2 n( , W H {l 1 z)oQ(z) 

^Q(z)=g , (11) 



we have set g = a/2l Q n^. 
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III. THE STOCHASTIC VOLTERRA EQUATION 



A. The Stochastic Volterra Equation and Its Solution 



The integral form of Eq. (jllj) is, 



Q{z) = Q z + g 



JO 



o Q(s ) ds ds 



On - // / ^-^-W H {l^s) o Q{s) ds. 



(12) 



Let us set the following initial conditions Q(0) = and Q(0) G Sjj. We are interested in 
finding a solution on the interval < z < L. What we have here is a stochastic Volterra 
equation with (Fredholm) kernel 



k H (z,s) :=g- 



which is continous and 



\k H (z,s)\\ H ^ q < gM X[o,z](s)s 1 \z - s\ 



(13) 



Now we have to find the conditions that make Eq. (fT2"|) solvable. We propose as ansatz 
the usual resolvent for convoluted kernels, that is, 

oo 

K H (z,s) = Y,K%\z,s), 



n=l 



such that 



Q{z) = QqZ + J K H (z, s) o (Q s) ds 



Qoo 



with the Kg given inductively by 



Kjg+%,8) 



K§\z,8)=k u (z }8 ). 



Kff (z,u)ok H (u,s)du, with n > 1, 



z+ I K H (z,s)sds 
o 



(14) 



(15) 
(16) 



0ksendaUi found this is the unique solution for bounded kernels in the distribution Hida 
space. The same procedure can also be applied to fractional Hida spaces. But it can be seen 
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from Eq. (jl3j) that it is not our case when s — > 0. Nevertheless, a new norm can be defined 
so that under it the former resolvent exists. Let it be, 



lit 



\LP>P'(J),-q 



sup / / \\G(z)k H (z, s)F(s)\\ H - q dzds 

pen -n<l J J J J 



\\LP(J),-q^ 1 - J J J J 



where ||i 7 '||LP(j),- g = (Jj \\F(s)\\' p H _ q dsY^ P is the Lebesque integral, and J = (0, L]. This 
norm is simplified using the Holder inequality the following relation can be proved: 

W 

•••^ii//.-,,' ''' 



\ kH \\LP,p'(J),- q < lmu 



J \JJ 



l^ljk H (z,s)\\ p H ^ q dzJ /P ds 



1/p' 



(17) 



Gripenberg et alQ discuss the deterministic counterpart of this construction. They proved 
that a resolvent solution exists whenever the norm of the kernel is less than one. This theorem 
can be tracked back to our norm in the stochastic case. Hence, the same hypothesis applies 



for this stochastic Fredholm kernel: \\k H \ 



LP'P'(J)-q 



< 1 for some q > 0. Then applying 



Eq. (JT3J) to the bounding condition (JT7J) we find 



\k H \\ L ^ iJ)} _ q < gMp- 1/p (%^V 
K h q \ sin 717/ 



< 1, 



since M is a small constant and g < 1. This guarantees the convergence of the proposed 
ansatz. 

The solution represented as a series of convoluted kernels, Eqs. ffTH) - (fTBj) . is useless for 
calculations. Next, we will prove that a fractional chaos expansion exists for the solution. 
Let us take the second term in the Wick product of Eq. ()14j) . it can be written 



X(z) =z + 



(n), 



Z, S 



n=l 



00 r pz 

s ds = z + / Kff\z, s)s ds 

n=l ^0 
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because it converges absolutely. The general term in this series can be written, 
K^\z, s)s ds 



1-H 




'Si 



Si 



oW H ( Sl )d Sl 




'0 Jsi 

; 1-// -.2 f Z f f Z T A n - 2 ), 



Si 




(s 3 - s 2 ){s 2 - si)(si - 0) 



J s\J S2 

z fZ l>Z 



S2S1 



52 51 



{gib'")*/ I I K%-' A} {*,**) — """" 1 "' ^'Wi'jIB'J 

\-LI ,i<l I I \ Z ~ S n)(s n - S n _i) • ■ • (Si - 0) H ^ ^ ^ H 



1-H\n 



Sri — 1 " 

n 



2 - S n 



n 



S-nSn— 1 ' °1 



"X[si_l,z) ( s i 



dB H ■ ■ ■ dB H 



«1 



igz H ll- H ) n f /W( Sn , . . . , ai ) dfl£ • • • dBj. 

</R™ 



We have used property (jHJ) to build the above adimensional integrals, and defined 

{Si — Sj_x) 



i=l 



(18) 



with s = 0. Now, the symmetrized form f^ n \s n , . . . , si) = (1/n!) Xlo-en /^( s o-„> ■ ■ ■ 5 s o-J 
induces the following relation 

/ . . . , ai ) dB% ■■■dB»=! fW(s n , ..., Sl ) dBZ ■ ■ ■ dB«, 

</r™ Jr™ 

and finally, 

x{z) = z\i + Y,( \rf (n \s n ,...,si) 

{ n=l M 1 

where g = log(z/l ) H . This will be nothing else but the fractional chaos expansion provided 

00 

Ef n i; ( ^<°° 



(19) 



n=l 



holds. In fact this condition express nothing else that the existence of the variance of the 
process, 



mi) = J 



i+$> 2 i/ { 



n)\2 



71=1 
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we used property (|B21[) . The search for an upper bound for the succession of 0-norms, given 
that the are symmetric, is straightforward: 

\f {n) \l = [ f [n \sn, • • • , ai)/ (B V„, • • • , ai)0(a„, <) • • • s[) ds n --- ds t ds' n - ■ ■ ds[ < 

< / / •••/ / (j){s ni s' n )---ct){s 1 ,s' l )ds n ---ds l ds' n ---ds' ll (20) 

JO J0 Js n -i Js' n _ 1 

because of definition (fTK|) and the fact < Si — < Si (idem < s' { — s^_ 1 < s^) the last 
inequality follows. Observing that 



1 r i 



\S n ., s n ^)ds n ds n 



n—1 J s n-l 



H(2H-1) / / \s n -s' n \ 2H 2 ds n ds' n 

^ s n. — 1 J s 'n — 1 

i - s n ^r + (i - s' n _ x f H 



< i, 



we iteratively apply it in Eq. (|2*U|) to find: |/' n ^|| < 1- Thus, the chaos expansion exists for 
all z < L whenever 

is satisfied. From the definition of g and the magnitude of the quantities involved in it we 
have: 

L«/;- 1/3H (C e 2 )- 1 / 2 ^. (21) 



So, the condition above is always fulfilled. 



B. Ray-light Statistics: An Example 

In this section we will use the stochastic ray-equation solution to study the statistical 
properties of the turbulent refractive index. Both coordinates of displacement are indepen- 
dent, and they also hold the same (non-coupled) differential equation. There is enough to 
consider the 1-dimensional then. The parameter election, Eq. (jA6|) . we have used in our 
treatment, also defines the meaning of the transversal velocities, for they are the angles of 
deviation. Being the velocities continuous we can set, 

Qo := IimQ(e) = 0| e=o G S* H . 
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Since our solution is dependent of the initial refractive angle 9, its behavior at the boundary, 
e — ► 0, must be known. This boundary is just the interface between turbulent and resting 
air. Henceforth, we will also model the initial angle as a fractional Brownian motion, 



9(e) = c[ X [oM s ) dBf = cB H {e), 



the constant c is adimensional and measures the strength of the noise. The length e works 
as a kind of correlation distance, as it goes to zero we are examining the properties of the 
interface's short-range correlation. 

Besides, any stochastic process can be put in terms of the spans described in the Appendix 
iBl and these depend on the construction of stochastic integrals by step functions. So, 
even if the former model needs to be corrected — maybe the interface introduces long-range 
correlations — the next results are useful; since, they are the building blocks for more complex 
stochastic processes. 

The solution (JUj) is written using the chaos expansion (fT^j) and the initial conditions: 
Q(z) = 9(e) o X z = zcB H (e) o ( 1 + f" ~g n [ f {n) (s n , ..., Sl ) dB^ ■ ■ ■ dB» ) . 

V n=l Jv + J 



From the Wick product properties is easy to see that 



E[Q(z)} = zcE[B H (e)] E 



oo „ 

i+vr / / (n) i 

n=l + 



s n ,..., Sl )dB*---dB£ 



-E[l] = 0. 



The evaluation of the variance from experimental data is the most common topic in 
many works related to the optical properties of turbulence because it is directly related to 
the Structure Constant. Hence, we calculate it using property (jB18|) . 



E[Q 2 



c 2 E(B»oX z y = c 2 E(D* X z ) 2 + E(X*)\ Xm 



E(X Z ) was already evaluated in the last section. The fractional Malliavin derivative appear- 
ing at the right-hand side demands elaboration, property (jB16|) implies 



Dth X 7 



Since the (^-differential is linear so we have 



(DfX z ) X[ o,e)( s ) ds - 



DfX z = z ~9 n D\ 



n=l 



f (n \s n ,..., Sl )dB"---dB 



(22) 
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We are going to compute these derivatives now: let us fix n > 2, from the first theorem {Bj) 
we can commute the stochastic integral and ^-differential, 



f^{s n ,...,s l )dB? n --.dBf l 



H 



= [ Df [ /» dBl ■ ■ ■ dBl dBl + / /<»> dBf n ■ ■ ■ dB»4>(s, Sl )d Sl . 

JR+ Jm™" 1 J JR™ 

Now, we recursively commute the operators, the 0-differential and the Wick integral. Each 
time we do so another integral as the last one on the right-hand side of the equation above 
is added. After {n — 1) iterations we reach the innermost integral, thus we evaluate 

/ f^\s n ,...,s l )dBf} = f f^(s n ,...,s 1 ) ( j ) (s n ,s)ds n , 
with the aid of property (jB17|) . Finally, 



[ f( n \s n ,...,s 1 )dB? n ...dB» =[ /» cf>(s,s n )ds n dB* n _ l .--dB? 1 + 
7r™ J Jr™ 

+ ...+ [ /» d Bl ■ ■ ■ fa, s k ) ds k ---dB* + ---+ [ /» dBf n ■ ■ ■ dB%<p(s, Sl )d Sl 
Jwt Jr™ 

/» 4>{s,s n )ds n 



n-l 
+ 



dB^-'-dB", 



to arrive to the last equality the symmetry of was employed. Instead, for n = 1 we just 
use property (jB17l) : 

/«(ai)d2J*| = / fM(s 1 ) ( f>(s,s 1 )ds 1 . 



Df 



Afterwards, we can build the fractional Malliavin derivative ()B15|) from the series 
D» X B =zg\ [ f^(s') Xm ( s )^ s ')ds'ds + 

oo „ „ 

+ V (n + 1) W / f {n+1) (s', ..., Sl ) X [o , e) (s)0(s, s') ds'ds 

n=l JV{.\J*% 



dB? ■ ■ ■ dBl \ ; (23) 



si [ ' 



its second moment is 
E (£>* X z ) 



2 ~2 

z g 



/ ( V) ds'ds 



+ 



+ in + lf~g 2n 



n=l 



f in+1 \s',-) X[0,e)( S )<t>(s,s') ds'ds 
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This series converges, we apply the same procedure as before to find a bound for the integrals. 
What is more, each norm appearing in the series is bounded by the zero term, 

2 



f in+1 \s',-)x{0,e)( S )^ S ') ds ' ds 

1 



{2(2if+l) ^ 



< 



.2H+1 



JO 



1 — s')<p(s, s') ds' ds 



-fl-e 



x2H+l 



2(2# + 1) L ~ " " J 2H + 1 " 2 I _ I' 

Thus, the existence of Eq. ()23|) is guaranteed. Finally, the variance of the displacements, 
with the norm 



|X[o, e )|^ = H(2H-1) I I \u - s 

'0 Jo 




2H-2 



duds = e 



2H 



is written as, 



E[Q 2 (z)} =z 2 c 2 



e 2H + f 



/ (1) ( s ') X[o,e)(s)(f)(s,s') ds'ds 



+ 



n=l 



~2n 



e 2H \f { % + (n + ir 9 2 



f (n+1) (s',-) X [O,e)OO0M') ds'ds 



Now, as the correlation distance goes to zero we recover the initial condition. While terms 
coming from the second moment of X z banish (they are all bounded and multiplied by e 2H ), 
it is not the case with those coming from the fractional derivative. We will not go through 
copious calculations since we are interested in a general outline of the solution; thereof, the 
solution can be expressed as 



E[Q 2 



Z C 9 . 2 2-2 ~2n/ 1 -t\ 

4 + z c g 2_^g (n + l; 

n=l 



We can estimate a bound for the second term: 



6^0 



F(~g 2 )=Aj2~9 2n (n + l) 



n=l 



f (n+1 \s',-) X[ o,e)( s )<t>( s , s ') ds'ds 



<4j2~9 2n (n + l) 2 (~) 

n=l ^ ' 

1 00 



n=2 
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Now, because x/(l — x) 2 = J2^=i nx2n ^ i s 

1 oo 1 oo 



r 



9 



n 2 ~g 2n - 1 



I 1 d 

2gdg 



g 2 



9 2 ) 2 



7.2 \ 



- 1 



n=2 n=l 

and then F(g 2 ) < g 2 / (1 — g 2 ) 3 , whenever g < 1. Finally, replacing the values for g, we have 



E[Q 2 (z)} 



Az 2 



2H+2j2/3-2H 







l + F[Az™l z Q 



2Hj2/3-2H 



(24) 



Furthermore, for the range of validity given in the past sections, the contribution of the 
function F is less than 10~ 6 . Thus, the first contribution to Malliavin derivative of X com- 
pletely characterize the variance once the interface's properties are defined. So, determining 
the behavior of the interface is crucial for the present model. 



IV. REMARKS & CONCLUSIONS 



In this paper we provided a fractional Brownian motion model for the turbulent index of 
refraction; afterwards, we used this model to build a stochastic ray-equation — a stochastic 
Volterra equation. For which we have given a (unique) solution. Our analysis just covers 
the H > 1/2 case which is meant for non-isotropic or near diffusive range turbulence; that 
is, the temperature gradients are relevant. 

Besides the fact that these Hurst exponents are hardly found in field experiences. They 
can easily appear within the laboratory. In particular, the turbulence is assumed per se 
nearly isotropic and homogeneous among optics works even in situations where the turbu- 
lence is created but not controlled. The example shown here can be used to test these cases. 
Also, we pretend to give a glimpse of the H < 1/2 problem extrapolating these results. 

We observe that the solution (|2*4^1 strongly depends on the initial conditions; our election 
of the initial angle at the example is the natural choice given the behavior of scalar quantities 
in turbulence. Under this condition, we apply Eq. (J24*j) to estimate the centroid's variance 
of a laser beam at a distance L — for any legth L less than 10 6 m according to Eq. ()21jl : 

V a r[Q(L)]^ C ^Air~ 2H L 2 + 2H , (25) 

where the correction to this result is of order O(F) ~ 10~ 5 . We must note that this is the 
contribution from the Malliavin derivative. That is, the constant term, of order zero, does 
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not contribute to the variance. Also, the anisotropy introduced by the termal flux should 
be observed in different values of c for each axis. 

Now, as H + — > 1/2 the variance of the displacement approaches to C^l^^L 3 . We notice 
that this behavior is the same found in Consortini et al, 1 . But it not comes from the 
isotropic Kolmogorov's law, with the hypotheses made at the introduction it corresponds to 
a Brownian motion (H = 1/2). That is, given a gaussian process with Structure Function's 
power as in Eq. (0) the result from Consortini does not follow. We have shown this can only 
be achieved on inertial-diffusive conditions, for ( n \ 1/3; shortly, the Beckman's hypotheses 4 
do not apply to gaussian models used in Fluids Dynamics. 

Nevertheless, the approach we have introduced is just the beginning of a long journey, 
since we must examine the stochastic ray-equation when H G [1/3, 1/2). By doing so, we will 
confirm that the power law for the variance of the displacements is univocally determined 
and no discontinuities arise, Eq. (|2*5j) is still valid. But, this will require the introduction 
of other tools — even new for the Stochastic Analysis — since in this range any smoothness 
property of the fBm processes is 1os1j^, and this inquiry will be the topic of future works. 

A step further should be considered afterwards. The proposed model for the turbulent 
index of refraction, Eq. (JJJ), just approximates the Structure Function. After we give a 
solution to the stochastic ray-equation for the whole range of Hurst parameters with this 
model, we will start examining other functionals of the fractional Brownian motion, which 
refines our theoretical results against the observed properties. 
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APPENDIX A: The Ray-Path Equations 



Fermat's Extremal Principle in Geometric Optics, means to find the variational solution 

to 

S[ I nds) = 0. 



The solution to this equation is interpreted as the ray trajectory, which we shall denote 
here as q{r). The parameter r has, in principle, no physical meaning. In Optics Treatises 
this parameter is usually replaced with one of the trajectory coordinates, which fulfills 
dqi/dr > 0, and is thus called the propagation direction. But the election of this parameter 
can not be done at will 1 -, since, for any parameterization chosen, the Optical Lagrangian 

L(q,q) = n(q)\\q\\ : 

(q, q G K 3 are the position and velocity respectively) is degenerated. Its solution is not 



univocally determined because - 



d 2 L 
dqtdqi 



for any pair (q,q). That is, calculating the momentum, 



Pi = M = n ^ q) uy (A1) 

we can write the Lagrangian as follows, 

r)T 

L (^) = E^=E^- ( A2 ) 

i i 

This is homogeneous in the velocities which allows us to recalculate the momentum and 
find, 



« =E f^ + y then 



so this matrix is singular as we stated above. Another consequence can be derived from 
Eq. (|A1|) ; also, it induces the following relation 

\\p\\ 2 = n 2 (q), 
18 



which indicates that the choice of coordinates and momenta is not free. 

Now, the degeneracy of the Lagrangian can be worked out in the Hamiltonian framework 
using the constraint we have just found. This problem of constrained Hamiltonians is known 
as Dirac's problem in the literature^. The procedure is as follows: first, define the constraint 
function 

*(p,?) = g Obll 2 - n2 (s)] ; 

then, calculate the Hamiltonian H from the original Lagrangian, from Eqs. (jAlj) and ()A2|) 
it is exactly zero 

H = J2 Pi ^ ~ L = J2 Pi ^ " YsP^ = 0; 

i i i 

finally, build a new Hamiltonian 

H(p, q) := H(p, q) + A* (p, q) = A* (p, q) 

and apply the variational procedure to all the coordinates included A. 
By doing so, we obtain the following dynamic equations 



. _ dH _ x &$ _ A v ^ 2 
^ dq dq 2 qU 
dH ay x 
op op 



(A3) 



and the constraint, over the phase space, 

= y(p,q) = ±[\\pf-n 2 (q)}. (A4) 

To ensure A is well defined we need to introduce compatibility conditions. Which arise 
from establishing the constraints as motion constants, that is, / = = {/, H} with {•, •} 
the Poisson bracket. In our problem these conditions reduce just to one: {H, ty} = 
which is automatically accomplished. There are no secondary constraints derived from the 
compatibility conditions so Eqs. (|A3|) and (|A4|) completely define our problem^. Notice 
that A is actually a smooth function on the constrained space that can be freely chosen. 
Finally, this pair of equations can be combined to yield 



(A5) 
I 2 = AV(g). 
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We now realize that with each selection we make for A the parameter r is also set, i.e. if we 
choose A = n~ l then 

||g|| 2 = 1 and ds = \\q\\dr = dr (A6) 

t is then the arc-length. But selecting A = 1 gives us ds = ndr and now the parameter is 
r = J ds/n. 

APPENDIX B: Fractional Stochastic Calculus 

Here we will introduce briefly the elements needed to build a stochastic calculus for B . 
The reader can find a complete reference in Hu and 0ksendal^ and Duncan et aZA 
Let H G (1/2, 1) be a fixed constant, and let us define 

<f)(s,z) = H(2H- 1) \s-z\ 2H ~ 2 , s,zGR. 

Then / G L^(R), if it is measurable and 

f(s)f(z)(f)(s,z) dsdz < oo. (Bl) 

Also the inner product can be defined in L|(R); 

(f,9h-= f(s)g(z)(f)(s,z) dsdz, for all /, g G Ir}(K), 

Jr Jr 

and it becomes a separable Hilbert space. 

Let 5(R) C Lj(R) be the Schwartz space of rapidly decreasing smooth functions on 
R. From its dual the probability space Q = <S'(R) can build, it is the space of tempered 
distributions to on R. Its associated probability measure, /i^,, can be found applying the 
Bochner-Minlos theorem, 



,2 




n 



where (u, f) is the usual pairing between elements in the dual and functions on R. With 
this definition it is easy to prove that 

E[(,/)] = 0, and E[(-,/) 2 ] = \f\\. (B2) 

The triplet (ft, B(Q) , n^) is thus a probability space — B(Q) is the Borel algebra on VL— 
usually called fractional white noise probability space. Then let L 2 ^^) = L 2 (Q,B(fl) 
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be the space of all the random variables X : Q — > R such that 

:=E|X| 2 <oo. 

Hence, those functions in L|(R) define the set of random variables of the form f(uj) = (to, f) 
which is included in L 2 (/^); that is, the condition ()B1|) induces square measurable random 
variables because of Eqs. (jB2j). 

It can also be shown that iS(R) is dense in L 2 (R). This means that for any / G L"^(M.) 
we can build a series f n 6 R such that f n — > f in £ 2 (R). What is more, the following limit 
exists in L 2 (/^): 

lim(u;,/ n > := (cuj). (B3) 

n— »oo 

We can now define the fractional Brownian motion process as follows: 

B H (z) := B H (z,u) = (u,X[o,z)) <= L 2 {^). (B4) 

In this definition B H is thought to be the z-continous version of the rightmost hand side 
term. The step function X\o,z) '■ ^ ~~ > [ — 1? 1] is defined by 



X[o,t]\ 



1 if < s < t 
-1 if t < s < 0. 
otherwise 



So taking / 6 L 2 (R), approximating it by step functions, and then using property (JB3 
traduces definition (IB4D to 



(w,/)= / f z alB H (z,u). (B5) 
It can be verified using the same procedure for /, g e L|(R) that 

E[(u;, /)<«;,</)] = (/,</),. (B6) 

Again let / be as above, as it is shown in Duncan et al. once defined 

£{f) = e W Uj dB H , (B7) 

the linear span 

S= \ j2a k £(f k );neN,a k eM,f k eLl(R) for k e {!,..., n} 



. k=i 
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is dense in L 2 (/x 9 i). 

Also there is another functional expansion from where this Lebesgue space with measure 
//^ can be build. This expansion is useful, in particular, to introduce some tools we will use 
through this work. To do so let us define the Hermite functions^ 

Ux) = {T-\n-l)\^y l/2 e- x2/2 H n ^(x)- n = l,2,--- 

on EL This set of functions is an orthonormal basis in L 2 (R). We can map this basis to an 
orthonormal one £ n = TT 1 ^ e L?(M), by means of the isometry 



rv(s) = c„ f (z-s) H -^ 2 f(z)dz, 

J \s,oo) 



[s,oo) 

where i 

H{2H-1)Y(\-H) 



It can also be shown that 



t(h - \) r(2 - 2H) ' 



/ | n (s) <j)(s, z) ds = c H I (z- s) H 3/2 £ n {s) ds, 
Jm J(-oo,z] 

because the £ n 's are an orthonormal basis these integrals are also smooth. 

Let X be the set of all finite multi- indices a = (a\, • • • , a m ) of nonnegative integers, we 
define 

H a (uj) := F ai ((o;,|i)) • ••H am ((w,£ m )). 

In particular, if we let e l := (0, • ■ • , 0, 1, 0, • • • , 0) denote the z'th unit vector then we get 
from Eq. (jB5|) and the definition of Hermite polynomials 

H e i(uj) = HHu,^)) = = [ 1(b) dBf. 

Jr 

These functionals are elements of L 2 (fi^), and they form its basis. That is, for X G L 2 (yU^) 
there are c a £l and a 6 X, such that 

X(cu) = J2c a n a (uj), (B8) 

and also 

W x \\h(^) = Y, alc ^ 
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These coefficients are given by c a = E[IH a ] ja\. With this property at hand we now define 
the fractional Hida spaces: the fractional Hida test function space Sh as the set of all 



ift(u>) = a a 7i a {uo) G L 2 (/i^,), such that 

« ! a K 2 ^) ka < °°, for al1 k e N > 



2 

H,k 



aeX 



where 



(2 N) 7 = JJ(2jp for any element 7 = (71, • • • , 7m ) G Z. 



The fractional Hida distribution space S H is the set of all formal expansions 

Y(uj) = '^^bpHpi^uj), such that 

pel 

\\Y\\ 2 H _ q = ^/?!aJ(2N) _9/3 < 00, for some q G N. 

pel 



(B9) 



With these definitions is not hard to see that Sh C £ 2 (/^) C S H . 

It is now time to show how the fractional white noise and integration with respect to B H 
can be defined. Let us first calculate the expansion for the stochastic integral in Eq. (|B5|) . 
For any / G L^(M.) — any given deterministic function — we have from Eqs. (jB8|) and ()B6|) : 

„ 00 

When / = X[o,z) i n the left hand side we recover Eq. (|B4|) and the following relation holds 



B H {z) =Y\\f ( [ £ k {s)<f>{s,u)ds) du 

k=1 U[0,z) \Jr J 



HA") e^, 



We can check its norm 



(2k)' q < 



l^)HSr.- s = E f / (f ik(s)<P(s,u)ds) du 
k=1 U[o,z) \Jr J 

< M 2 z 2 2-o jr k 1 '^ = M 2 z 2 2-% (q - , 



(BIO) 



(( is the Riemann's zeta function) this because of 



i k (s)(f)(s,u) ds 



[u — s 



,H-3/2 



€k(s) ds 



(— oo,u] 



< Mk 1/6 



23 



according to Szeg When q > 4/3 the former inequality shows that B H is continuos and 



differentiable in S^. Hence, 



(uj) ■= w H (z) e S 



H 



(Bll) 



is the formal definition of fractional white noise. Again it is also continous, when z > s 



\W H (z)-W H (s^ 2 



H-q 



5>' ! 

k=l 

00 



£ k (u)(f)(z,u)du- / £ k (u)<f>(s,u)du 



{2k)- q < 



k=l 



n 2 



[(^_ s )_ w ]^-3/2|^( s+u )| du 



(2k)' q < 



<2-ic 2 H M 2 ([q + - 



[{z -s)-u\ 



H-3/2 



2WY l\ ,,, 



(B12) 



We now define the Wick product, like in Hu and 0ksendal, as follows: let be X(uj) 

EaeJ a «^«(^) alld Y ( U ) = E/JeX in S H then 



a,/3SX 7GX \a+/3=7 / 



(B13) 



this product is commutative, associative and distributive like the usual product in R. Under 
the norm || ■ \\n-q it can be shown that all the power series defined in the real line have 
their counterpart in the distribution Hida space. If we understand by X on the n-times Wick 
product of X 6 S^, then the Wick exponential is 



exp 



00 1 

w = E ~ x ° n - 



n=0 



In particular when X = (u, f) = J R f s dB^ holds, then 

exp*«w, /» = £(/). (B14) 

The right-hand side was defined in Eq. ()B7|) . and therefore we have shown the link between 
the two representations presented up to now. 

It is appropriate to introduce the fractional Malliavin derivative or 0-derivative, for X 6 
L 2 (^) and g G Lj(R), 

D $9 I(w) = limMx(cj + 5 A ($g)(u)du))-X(u) 
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where ($g)(z) = J R 4>(z,u)g u du. Afterwards, if there exists a function DfX such that 

D„ g X = [ (DfX)g s ds, Wg G Lj(R), (B15) 

we say that X is 0-differentiable, and DfX is the ^-differential. These are differential 
operators, and they also present the following properties: let X be as always and /, g : R — > R 



then 



D* B (u,f) = (f,g) 4 ,, (B16) 
D+(u>,f)= [ (f>(u,s)f u du (B17) 



Moreover, it can be proved for any two 1,7 G L 2 (^) — decomposed by the span £ — using 
the Wick product and the 0-derivative properties that, 



E 



Xo f f s dB? Yo [ g s dBf 



E[(D^ f X)(D^Y)+XY(f,g)^\ . (B18) 



These properties allows to change the integrator inside Eq. (jB5|) by a stochastic function 
X : R x Q — > R. The basic procedure consists of building a Riemann sum, replacing the 
standart product by the Wick one, 

n-l 

S n (X) := ^X 2i « +i - B%), (B19) 

i=0 

and finally observing its convergence (whenever E |X|^ < oo) under the E |-|^-norm. More- 
over, the following equality holds 

X s dBf= I X s oW s H ds, 

[Q,L) J[0,L) 

while the integral on the left-hand side represents the limit of Eq. (|B19|) . the right-hand 
side is just the integral evaluated under the Hida expansion of the Wick product defined in 
Eqs. ip5 j) -l|BI5 ) . 

We will complete this Appendix enumerating some useful theorems from Duncan et al. 27 : 
Let X z a stochastic process defined as above, and sup^ g [ L ^ E |D^X|^ < oo. Also, let 

= l[0,z) f» dB s ■ Then f° r S ' Z e [°> L ) 



Df Vz = [ DfX u dB^ + / X u <f)(s, u) du. 

J[0,z) J[0,z) 
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Let r\ z = J, z) f s dB^ and F(z, x) : M.+ xl-tl, where f is continous and F has second 



l[0,z) 

continous derivatives. Then 



F(z,r ]z ) = F(0,0)+ / —(8,<o.)d8+ / —(s, Vz )f s dB? 

J[0,z) ° s J[0,z) OX 

f d 2 F f 
+ / "aT^Vs) / <j>(s,s')f s/ ds'ds. 

J\0.z) ox J\0.s) 



[0,z) ux J[0,s) 

2 



If X G L 2 ^^) then there exists a sequence {f n G L^(R")} nG N such that Yl™=i 1/nL < 00 
and 

00 „ 

X = E[X] + J2 fn(su • ■ ■ , s n ) dB» ■ ■ ■ dB? n . 

n=l 

Here it is understood L|(R") as the n-dimensional space of symmetric functions and 
\fn\l = / /n(si, • • • , s n )f n (s[, • • • , s' n ) <f>(si, s[) ■ ■ • 0(s n , s' n ) ds x ■ ■ ■ ds n ds[ ■ ■ ■ ds' n . 

L|(M") is the n-dimensional space of symmetric functions. Given the base complete or- 
thonormal base {£ n }„ e N C L^,(R + ) then L|(R") is the completion of all function of the 
following form: 

f[si,...,s n )= ^ a k 1 ,...,k n iki(si)ik 2 (s 2 ) ■ ■ -^ kn ( s n)- 
i<fci,...,fc n <fc 

Associated to this functions we define the multiple integral as 

l<kl,...,k n <k J R ^ R ^ M 

It is not difficult to prove 27 (Lemma 6.6) that given / G L|(R") and / G Lj(R^) it is 

{(/, o"L if n = m 
(B21) 
if n 7^ m 

Moreover, for the iterated integral 

J0<Sl<S2<-<S„<t 

/* f / fn(si, " " " , Sn-1, «n) dB% • • • dfljO dfl£ 

JO \V0<Sl<S2<-<S n _l<S„ / 

is n! times /„(/). 
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